function M = madist(mdl, X, op)
%MADIST Computes the Mahalanobis distances to the Gaussian center
%
% [ Syntax ]
%   - M = madist(mdl, X)
%   - M = madist(mdl, X, 'shared')
%
% [ Arguments ]
%   - mdl:      the model or an array of models
%   - X:        the sample matrix
%   - M:        the resulting matrix of distances
%
% [ Description ]
%   - M = madist(mdl, X) computes the distances from the samples in X
%     to the center of Gaussian model(s).
%
%     Suppose there are n samples, then X should be a matrix of size
%     d x n, with each column giving a sample. 
%
%     If mdl is a single model object, then in the output, M is a row
%     vector of size 1 x n. If mdl is an array of K models, then M is
%     a row vector of size K x n. 
%
%   - M = madist(mdl, X, 'shared'). If you know that all models in mdl
%     share the same covariance matrix, you can use this syntax to tell
%     the function, such that it can make use of this fact to accelerate
%     the computation.
%
% [ Remarks ]
%   - If mdl is an array of multiple models, then all models should be 
%     in the same cache level. 
%
% [ History ]
%   - Created by Dahua Lin, on Dec 25, 2007
%

%% parse and verify input

assert(isnumeric(X) && ndims(X) == 2, ...
    'sltoolbox:slcisogauss:madist:invalidarg', ...
    'X should be a numeric matrix');

K = numel(mdl);
d = getdim(mdl(1));

assert(size(X,1) == d, ...
    'sltoolbox:slcisogauss:madist:invalidarg', ...
    'The dimension of samples does not match that of the models.');


if nargin >= 3
    assert(ischar(op) && strcmp(op, 'shared'), ...
        'sltoolbox:slcisogauss:madist:invalidarg', ...
        'the 3rd argument should be a string "shared"');
    
    shared = true;
else
    shared = false;
end




%% main

mu = [mdl.mu];

% terms

txx = sum(X .* X, 1);

if K == 1
    tuu = mu' * mu;
else
    tuu = sum(mu .* mu, 1)';
end

txu = mu' * X;

% combine

if K == 1
    M = txx + tuu - 2 * txu;
else
    M = bsxfun(@plus, bsxfun(@plus, (-2) * txu, txx), tuu);
end

% scale

if K == 1 || shared
    s = 1 / (mdl(1).sigma);
    if s ~= 1
        M = M * s;
    end
    
else    
    s = 1 ./ ([mdl.sigma])';
    if ~all(s == 1)
        M = bsxfun(@times, s, M);
    end
    
end

        









